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PREFACE: 


A vectorized Poisson solver, specifically designed 
for the CDC 200 computer series, was developed to obtain 
solutions of self-adjoint differential equations of the 
type : 


a(p) V • [b(<p,A,p) Vu] + c(<p,\) 3 2 pU = source (<p , A , p ) 

< V : = V <pA > 

where u is the unknown variable and a, b, c represent 
the coefficients as functions of the spatial variables <p , 
A, p . The purpose of this solver is to generate solutions 
for large gridsizes with spherical geometry at a high 
efficiency rate. It is assumed that the radial extent of 
the physical domain is much smaller than its average 
distance from the center of the sphere. 

The two -parameter Chebyshev method has desirable 
properties such as an effective convergence rate and simple 
implementation for vector processors. The algorithm can be 
applied to a variety of problems with spherical geometry. 

It allows to obtain solutions of the quasi-geo- 
strophic omega- equation if its coefficients are spatially 
distributed such that the equation remains elliptic. A 
three-month G CM simulation provides the data base for the 
computation of the quasi -geostrophic omega solution which 
is verfied against the fields obtained from the divergent 
wind fields . 

The results presented in this publication are part of 
a more extensive diagnostic study of the mechanisms 
responsible for the maintenance of the general circulation. 
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INTRODUCTION 


The development of a vectorized Poisson solver was 
motivated by the need to obtain solutions of the quasi- 
-geostrophic omega- equation. 

With the aid of the quas i - geo s t r oph ic omega - equat ion 
the driving forces of the general circulation can be 
examined individually since partial solutions of the 
omega- equation are obtainable as response to individual 
source terms. Therefore the relative roles of dynamical 
and physical processes and orography in producing the 
vertical velocity distribution can be ascertained. 

Primitive equation models provide a nonlinear omega 
from the mass continuity equation using the divergent wind 
fields. The "primitive" omega is more accurate, and to the 
extent that quas i - geo s trophic theory is valid, compares 
with the linear quas i - geo s trophic omega. The effort to 
calculate the quas i - geo s tr oph ic omega is much greater, and 
thus requires the development of an efficient three-dimen- 
sional Poisson solver which is applied to a spherical 
grid of 72 x 46 x 9 gridpoints. 

In Section 2, we express the analytic form of the 
omega- equat ion in flux form to facilitate the application 
of a suitable finite difference scheme. The physical 
domain is geometrically a spherical shell with large 
radius. In Section 3, the finite difference equations 
are developed, using the box method which preserves 
symmetry of the finite difference equations. Section 4 
describes the vector layout of a point Chebyshev solver for 
CDC 200 series computers. The convergence rate is 

accelerated by the choice of two iteration parameters 
which are functions of the spectral radius of the iteration 
matrix. In Section 5, we present a discussion of some 
results from a three-month integration of the GLA 
Fourth-Order General Circulation Model (see Kalnay et al . , 
1983). It includes a comparison of the "primitive" model 
omega field and the recalculated quas i - geos trophic omega 
field, which are verfied against the observed "primitive" 
omega from FGGE III b data. 


FORMULATION OF THE QUAS I - GEOSTROPHIC OMEGA-EQUATION 
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The vorticity, divergence, thermodynamic, and 
continuity equations describing the linear balanced model 
are : 


d t V 2 V> + • V ( £ + f ) + • Vf = fdpW (2.1) 
fV 2 ^ + Vf • VV> - V 2 $ (2.2) 
3 t 3p<£> + v • V3p$ + au> = -Q (2.3) 
V 2 * + 3 p w =0 (2.4) 


where Q denotes the sum of all heating terms and a 
represents the static stability parameter (see Haltiner and 
Williams , 1980) . 

Upon multiplying (2.1) by f and introducing (2.2) 
into (2.1), the vorticity equation becomes: 

3 t V 2 $ - Vf • V3 t ^ + fv^ • V(£ + f) + fv^, • Vf = f 2 3pW 

(2.5) 

Upon eliminating time derivatives between equations 
(2.3) and (2.5) [with the aid of (2.2)] one obtains the 
omega -equation: 


V 2 (cr«) + f 2 3 2 p w = f0 p [v^ • V(£+f)] - V 2 (v • V3 p $) 

- V 2 Q + fVf • 3 p v x - Vf • V3 p 3 t V> (2.6) 


It takes the form of the quas i - geos trophic omega- equation 
if one neglects the last term in the foregoing relation- 
ship. By virtue of scale analysis this term can be shown 
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to be negligible compared with the differential vorticity 
advection in midlatitudes and compared with the heating 
in the tropics . 


Without introducing any sensible error we may also add 
the small term f^p( v y * V£) which [with the neglect of 

the last term in (2.6)] enables us to write (2.6) in the 
form : 

V 2 (aco) + f 2 a 2 p u> = f3 p [v • V(£ + f)] - V 2 [v • V(3 p $+Q)] 

(2.7) 


which is desirable for reasons of computational simplicity. 

The parameters a and f slowly vary over the 
globe. Following Hoskins (1978) we may treat the static 
stability parameter as a constant over length scales of 
the order of the Rossby radius and smaller, but as varying 
over larger scales. To ensure the latter assumption, the 
static stability is computed as an equally- weighted nine 
point horizontal average in space. It is therefore a 
function of all three spatial variables whereas the 
Coriolis parameter naturally depends on <p only. As a 
result, the elliptic omega- equation can be approximated 
with little error in flux form: 


3^(a/9 r d \a>) + cos (p d ^(cosep a/0 T 

+ r 2 cos 2 <pd 

— f r cos<p 3p[u 3^(£+f)] + f r cos 2 <p 3 p [v 
+ a/e [d 2 A + coscp d (p (coscpd cp ) ] 

[u/(r cos<p) 3^0 + 

- a/0 [3 2 a + cos<j p dpicoscpd'p) ] Q 


p(f 2 3p«) 

v* +f > ] 


v/r a^e] 


( 2 . 8 ) 
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where T 


* -3p0 denotes the static stability. 

The equation can be rewritten concisely as: 

L w (w) := <(V , 3 p ) , (a/9 T V,f 2 3 p ) > w 

; - f 3 p (v • V($+f) + a/e V 2 (v • ve - Q*) 

(2.9) 

! where Q* - 9/a Q and < a , b > denotes the scalar product 

of a and b . In this form, it becomes apparent that the 
associated operator L w is self-adjoint in an integral 
sense with respect to the weighting function l/(a/e T f 2 ). 

The boundary conditions at the bottom and the top of 
the atmosphere are given by the expressions 

w “ -Ps v s *V($ s -$) 

(where for our climatological analysis v s • is neglect- 

1 ed in consistence with our quas i - ge o s tr oph ic assumption), 

and 

I 

i 

co * 0 , respectively. 
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3. FINITE DIFFERENCE FORM OF THE OMEGA- EQUATION 


Above formulation of the partial differential 
equation in flux form allows symmetric finite difference 
equations although the grid need not be uniform. This 
greatly reduces storage as shall be elucidated in this and 
the following section. 

Equation (2.9) of Section 2 is integrated at each 
mesh- point over a box with the meshpoint at its center: 


(v,a p ) • (a/e r v,f 2 a p ) u dv - 

V S (V) 


RHS dV 
V 

(3.1) 

where RHS represents the right-hand- side of (2.9) and the 
volume integral of the lef t - hand - s ide has been converted 
into a surface integral with S the surface of volume 
element V . 


E (rear) 


N 

dA 

rdp W 

Figure 1: geometric configuration 



(front ) 




(a/e r v w ,f 2 a p w).ds 
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The volume element is given by dV = r 2 cos cp dcp dA dp 
and the surface elements of the faces are: 

r cose p dA dp in the South-North (S-N) direction 

r 2 costp dA dcp in the bottom-top direction 

r dcp dp in the West-East (W-E) direction 

To obtain the finite difference equations, we introduce the 
indices i, j, and k denoting W-E, S-N, and bottom-top 
directions running from 1 to I , 1 to J , and 1 to K 

respectively. Thus, we obtain: 

» A A 

RHS dV - RHS c r 2 cos 2 <p c A A A <p Ap (3.2) 

V 


and the different face contributions to the coefficients 
from the surface integral are: 

upper face: f 2 ( w k+ 1 " / (Pk + 1 “Pk^ r2 cos2( P c A<p AA 

lower face: f 2 ( w k - 1 " w k) / (Pk~ Pk - 1 ) r2 cos 2( Pc AA 

front face: a front ( w i - 1 ‘ w i ) / ( r cos <p c AA) r A<p Ap 

rear face: a rear ( w i + 1 “ 00 i ) / ( r cos <p c AA) r A cp Ap 

Northern face: (coj +1 " w j ) / ( rA<p) r cosep^ AA Ap cos cp c 

Southern face: ag ( u> j _ ^ - u> j ) / ( r Acp ) r cosipg AA Ap costp c 

(3 . 3-3 . 8) 


where subscripts c, S, N indicate evaluation of the corres- 
ponding variable at the center, the Southern face, the 



Northern face of the volume element, respectively, and Ap 
is defined as : 


1/2 (Pk+l"Pk-l) for interior nodes. 

a front> a rear» a N anc * a R are obtained by linear 
interpolation from neighboring meshpoints . 

In finite difference form, the source terms and LHS 
omega terms are computed consistently with o r der ( di s t anc e ) 
truncation error. Periodic boundary conditions are speci- 
fied in the East-West direction, Dirichlet conditions at 
the bottom and top of the atmosphere which are assumed 
to be at lOOOmb and 50mb. 

The kinematic condition at the bottom expands to: 


w - -P S / R T S ( u /( r COS <p) a A $ s + v/r d^s) (3.9) 

The grid is treated as a rectangle with multiple 
representation of the North and South poles. Half -box 
geometry is applied along the six boundary planes. The 
coefficients across the boundary faces of these half-boxes 
are taken to be zero (zero-flow condition). 

The periodicity condition along the meridian boundary 
at 180° longitude East and West is achieved by adapting the 
following procedure: 


Half-box geometry factors are applied at both boundaries 
(180° West, 180° East) to obtain the coefficients. During 
each iteration the two separate corresponding solutions 
pertaining to the half-boxes are arithmetically combined to 
form a single solution (see algorithm in Section 4). 

The spatially symmetric formulation of the finite 
difference equation requires only storage of four coeffi- 
cient sets: a set of diagonal coefficients and one set for 

each geometric direction. These coefficients are specified 
on a three dimensional grid with index i starting at 180° 
West, index j starting at the South pole, and index k 
starting at the bottom of the atmosphere. 


7 


assigning a large number, say 10^, to the diagonal 
coefficient of those meshpoints on the boundary and outside 
the physical domain. Correspondingly, the RHS must be 
specified as the product of the unknown variable u> and 
this pre-set large number on the boundary. 
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4. DESIGN OF A VECTORIZED POISSON ALGORITHM FOR CDC 200 
SERIES 


It is important to select a suitable layout for 
indexing the arrays and a method which allows vec tor iz a t ion 
of the iteration scheme. 

As pointed out in the foregoing section, special 
consideration must be given to the boundary nodes which 
include the polar nodes, the nodes at the top and bottom of 
the grid and the meridian boundary where the solution 
from both ends must be joined: 

The adjacent longitudes outside the physical domain 
are occupied with a set of zero off-diagonal coefficients. 
This ensures vec tor izat ion of the entire code without the 
need for checking indices in any of the three geometric 
directions during the iteration process. 

Six descriptors "D" for the six different directions 
and one descriptor for the center of the individual meshbox 
pointing to the proper addresses of the coefficient vectors 
are introduced: 


D(CW) - D ( CE ) 

- 1 

D ( CS ) = D ( CN) 

- I 

D ( CB ) = D ( CA) 

- IJ 

D ( CP ) 



where CE/CW, CN/CS and CA/CB are three off-diagonal 
coefficient sets referring to the East, North, and top 
directions, CP the diagonal coefficient set and " D" 
descriptors pointing toward the proper locations. 

Similarly, six descriptors plus an additional one are 
also needed for the solution vector of the unknown variable 
0) 

As iteration procedure, a point Chebyshev method with 
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adaptive determination is implemented because of its ver- 
satility and applicability for solutions of large three- 
-dimensional elliptic equations. This iterative method can 
be ideally vectorized and is proved to be efficient for 
large 3D elliptic equations. 

Two iteration parameters, alpha (a) and beta (/3), 
need to be specified which are updated prior to each 
iteration step. The rate of convergence depends largely on 
the choice of these parameters. The fact that the number 
of vertical planes is much smaller than the number of 
planes in any of the two horizontal directions leads to 
rapid convergence. For the present problem, we select a 

spectral radius of 0.9 which detemines the coefficients of 
the Chebyshev recursion formulas for the iteration 
parameters : 


Or ( 1 ) 

:■= 0.95238 

a(2) : - 

1 .6123 



a(t) 

:= 1. / [1.05 

- 0.225625 

a(t-l) ] 

for 

t>2 

>3(1) 

: = 0. 
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M 

0 

U1 

Q 

/-N 

rt 

1 

1 . 


for 

t>l 


(4.1) 


With these choices, alpha approaches a value of approxi- 
mately 1.336 and beta approaches a value of approximately 
0.403 as t approaches infinity. 

Upon defining UNOW as the solution vector, UM1 and 
UM2 as auxiliary iteration vectors, SOURCE as the RHS and 
CPINV as the vector of reciprocal diagonal coefficients, we 
can summarize the iteration algorithm in CDC Vector Fortran 
as follows : 


UNOW * -SOURCE 
t = 1 
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100 UNOW ( i=l , j , k) = UNOW ( i=I , j , k) replaced by 

UNOW ( i =1 , j , k) + UNOW ( i=I , j , k ) 


(to account for periodicity condition) 
UNOW = CPINV*UNOW 

UNOW = ALPHA(t)*(UNOW - UM1) + BETA(t)*(UMl - UM2 ) 
USPOLE(k) = average over i of UNOW along line (j=l) 
for k=l , K 

UNPOLE(k) = average over i of UNOW along line (j=J) 
for k=l , K 

UM2 = UM1 

UM1 - UNOW + UM1 

UNOW = -SOURCE + CE*UM1E + CW*UM1W 

+ CN*UM1N + CS*UM1S 

+ C A*UMl A + CB*UM1B 

(where appended letters E, W. . . indicate addresses 
pertaining to the East, West. . . directions) 

t = t + 1 

go to 100 (until convergence is reached) 

UM1 and UNOW are periodically compared in order to 
terminate the iteration when sufficient accuracy is 
achieved . 

To achieve an accuracy of 10~ 5 starting with the 
highest significant digit, approximately 300 iterations and 
1-1.5 CPU seconds on a CDC 205 are needed for our appli- 
cations to a grid with approxiamtely 40,000 gridpoints 
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and given source distribution (The value 40,000 is 
ideally below the critical limit for vector processing on 
the CDC 205, namely 2 16 - 1 = 65,535). 

The vectorized algorithm is remarkably superior to 
its scalar, but otherwise identical analogue with a CPU 
performance of approximately 100 CPU seconds. 

The above described algorithm is applicable to a wide 
variety of other Poisson problems with spherical geometry 
which, for example, frequently arise in the context of 
electrostatics and heat flow. 
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5. EXAMPLE OF AN APPLICATION OF THE VECTORIZED POISSON 
AND DISCUSSION OF RESULTS 


5.1 Observed and GLAS-GCM omega fields 

Figure Al shows a three month average of the 500mb 
vertical velocity (omega) field for the period from Dec 15, 
1978 to Mar 15, 1979, calculated by Krishnamurti and Sheng 
(1984) from unitialized ECMWF FGGE Illb divergent wind 
field patterns. Bands of rising motion (unshaded areas) 
are clearly centered near the equator and sixty degrees 
North and South latitude with bands of sinking motion 
centered around thirty degrees North and South. An intense 
center of sinking motion is found on the southern slope of 
the Himalayas north of the Bay of Bengal. Another 
prominent feature is the center of rising motion extending 
south-eastward from South America almost connecting with 
the ascending belt around sixty degrees South. There is 
also an intricate vertical motion pattern with rising 
motion across the Andes. 

Figure A2 shows the corresponding vertical motion 
pattern derived from a GLAS three-month simulation using 
the formula 


m * 

u) an + an (5.1) 

where a is the vertical velocity in a coordinates, 
obtained from divergent wind field pattern, and n is the 
surface pressure. This field is interpolated from the 
model grid to a pressure coordinate system and then to the 
FGGE horizontal grid by means of bicubic splines. As in 
the case of the vertical motions derived from the ECMWF 
FGGE Illb analysis, the simulated vertical velocity pattern 
exhibits bands of sinking motion in the proximity of thirty 
degrees North and South. It also displays a band of rising 
motion extending SE from Brazil to the rising band at sixty 
degrees South. The smaller scale centers of vertical 
motion are, however, more intense than those of Figure 
Al and do not correlate well with the latter. The corre- 
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lation coefficient between the area-weighted omega fields 
at 500mb shown in Figures A1 and A2 is in fact only 
0.27 + 0.03, which means that the simulated field accounts 

for only about 9% of the variance of the analyzed field. 

A notable feature of the simulated field is the 
bimodal distribution of vertical motion associated with 
each of the major mountain chains, with East-West 
orientation in the case of the Rocky Mountains, the 
Himalayas and the Andes, and a North-South orientation 
in the case of the Greenland Plateau. (The GLAS enhanced 
topography field is shown in Figure A3.) The vertical 
motion pattern in the vicinity of the Andes is more 
complicated than in the FGGE analysis. There is also a 
prominent region of the upward motion in equatorial East 
Africa which has no counterpart in the FGGE analysis. 


5.2 Quasi-geostrophic omega field from GCM data 

In order to investigate the features of the GLAS 
simulated model omega we solved equation (2.9) with the aid 
of the vectorized Poisson solver to obtain the quasi- 
-geostrophic vertical velocity field. Scale analysis 

of the most general form of the ome ga - e qua t i on reveals that 
the quasi-geostrophic version is valid in middle and high 
latitudes where the Rossby number is small. The same 
analysis reveals that near the equator the dominant terms 
are VT • Vw and Q , with all other terms being negli- 
gible. Since the quasi-geostrophic omega equation also 
reduces to a balance between these two terms as the 
equator is approached, we shall use it rather than the more 
general omega- equation arising from the balance equation 
which is more cumbersome to solve. 

To the extent that the quasi-geostrophic omega- 
-equation is valid and friction forcing has a negligble 
effect on the vertical motion at 500mb, the solution of 
equation (2.9) forced by the sum of all other forcing terms 
on the right should resemble the omega field obtained from 
the GLAS model simulation. Figure A4 reveals that it does, 
although the maxima and minima are somewhat weaker than 
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those in the original simulation. Figure A4 was obtained 
by time - averaging the individual solutions of the omega- 
-equation at 181 twelve-hour intervals from 00Z Dec 15, 
1978 through 00Z Mar 15, 1979. This is necessary in 

principle because static stability, which appears as a 
coefficient on the LHS of the ome ga - e qua t i on , varies with 
time and could conceivably correlate with go . 
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8. APPENDIX 


The following figures referenced in the text show the 
observed FGGE omega field at 500mb, ^Figure Al), the GLAS - 
-GCM model omega field (Figure A2), the bottom topography 
(Figure A3), and the quas i - geo s trophic omega field (Figure 
A4 ) . 

The contour interval in Figures Al , A2 and A4 is: 

0.5 1 0 " ^ mb / s . 

The contour interval in Figure A3 is: 

0 . 5 10 4 m 2 / s 2 . 
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